Free energy calculations: An efficient adaptive biasing potential method 
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We develop an efficient sampling and free energy calculation technique within the adaptive biasing 
potential (ABP) framework. By mollifying the density of states we obtain an approximate free energy 
and an adaptive bias potential that is computed directly from the population along the coordinates 
of the free energy. Because of the mollifier, the bias potential is "nonlocal" and its gradient admits 
a simple analytic expression. A single observation of the reaction coordinate can thus be used to 
update the approximate free energy at every point within a neighborhood of the observation. This 
greatly reduces the equilibration time of the adaptive bias potential. This approximation introduces 
two parameters: strength of mollification and the zero of energy of the bias potential. While we 
observe that the approximate free energy is a very good estimate of the actual free energy for a large 
range of mollification strength, we demonstrate that the errors associated with the mollification may 
be removed via deconvolution. The zero of energy of the bias potential, which is easy to choose, 
influences the speed of convergence but not the limiting accuracy. This method is simple to apply to 
free energy or mean force computation in multiple dimensions and does not involve second derivatives 
of the reaction coordinates, matrix manipulations nor on-the-fly adaptation of parameters. For the 
alanine dipeptide test case, the new method is found to gain as much as a factor of ten in efficiency 
as compared to two common adaptive biasing force formulations and it is shown to be as efficient as 
well-tempered metadynamics with the post-process deconvolution giving a clear advantage to the 
mollified density of states method. 



I. INTRODUCTION 

Many interesting physical systems can be categorized 
as rare-event systems. The uniting feature of these sys- 
tems is that the dynamics involved require a time reso- 
lution much smaller than the timescale on which inter- 
esting events take place. Of central importance to the 
evolution of such systems is the free energy. Low ly- 
ing regions of the free energy and the barriers separating 
these regions dictate the thermodynamics and, to some 
extent, the kinetics Q of the system. Because free energy 
barriers are only rarely crossed, efficient exploration of 
the free energy landscape is practically impossible with 
straightforward integration of the equations of motion. 

Recently, a number of approaches have used the his- 
tory of the dynamics to accelerate exploration of the 
free energy landscape 0-0] • I n t nese methods informa- 
tion about the free energy is estimated during simula- 
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tion and that information is fed back to the dynam- 
ics as a statistical bias. While there are many varia- 
tions on this idea, the common aim is to minimize the 
time spent sampling regions of the free energy that have 
been sampled in the past. These schemes may be clas- 
sified into two categories: adaptive bias force methods 
(ABF)@, i] which use an approximation of the mean 
force to bias the dynamics; adaptive biasing potential 
methods (ABP) 0, 1,0 

which use an approximation of 
the free energy as a bias potential. 

The underlying idea for all adaptive methods is that 
it is computationally more efficient to sample the distri- 
bution associated with a flattened free energy than it is 
to sample the density associated with the actual, very 
rough free energy. We propose an ABP method that 
builds an approximate density of states (DOS) and uses 
that approximation to define a bias potential. Mollifi- 
cation of the underlying density of states produces the 
desired approximation and leads to a smooth, adaptive 
bias potential whose gradient admits a simple analytic 
expression and that can be computed without knowledge 
of the actual density of states. Because the actual and ap- 
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proximate free energies are related by a convolution, it is 
easy to recover the former from the latter via deconvolu- 
tion. Our framework is not restricted to one-dimensional 
or orthogonal reaction coordinates. Moreover, it avoids 
second derivatives of the reaction coordinate. 

This paper is organized as follows. We describe our 
ABP method in Section [XT] (see in particular Eqs. ©, 
(fTU)l and (dU), and comment on its convergence. We 
contrast this method to existing ones in Section ITO1 and 
present some numerical validation on a benchmark sys- 
tem in Section IIV1 Our conclusions are summarized in 
Section M 



II. DESCRIPTION OF THE METHOD 

A. Free energy and its mollified version 

Consider a system whose configuration is described by 
a variable x G X, where X is the configuration space. 
We denote by V the potential energy function. Assume 
that we are given a .ZV-dimensional reaction coordinate 
£(x) with values in f2, which characterizes some physical 
event. The density of states at a value £* of the reaction 
coordinate is defined as 

e -M«*) =z -i f 5{^{x)-C)e- &V{x) dx, (1) 
J x 

where ft = l/(fegT), and Z is a normalization constant, 
chosen such that 



-M(£) 



d£_ = 1. 



Eq. (TT|) defines the free energy A(£*). Recall that, in 
practice, the free energy needs only be known up to an 
additive constant since the important quantities to de- 
scribe the relative likelihoods of physical states are free- 
energy differences. 

In general, the free energy is unknown and has to be 
approximated. The method we propose in this work is 
based on the following limit: 



e-^«*> = hm e -^«(n 

a->0 

where 

e -M<»(n = z- 1 [ 6 a (£(x) - D e- fjV{x) dx, (2) 
J x 

with for instance a Gaussian approximation of the Dirac 
delta function: 
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Equation ([2]) defines an approximate free energy A a , ob- 
tained by sampling the density of states at a finite a, i.e., 
by sampling a mollified density of states. Notice that 



the approximation resulting from finite a can in fact be 
rewritten as a convolution of the actual density of states 
e~@ A with S a . Indeed, 



e -Mo,(n 

= Z~ l f 5 a (ax)-C)e- mx) dx 
J x 

= Z~ l [ [ *a(f-0 -IK 
Jn J x 



(3) 



-f3V(x) 



dx d^ 



= / s a #-e)e 

Jn 
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This remark is the basis for an extraction of the actual 
free energy A from A a through a deconvolution proce- 
dure (see Section IIV A[) . While we make this presenta- 
tion with a scalar a, this could easily be generalized to 
the case where a takes different values in different dimen- 
sions of the reaction coordinate. 

Equation © is also helpful in assessing the errors in- 
troduced in umbrella sampling (US) and thermodynamic 
integration (TI) simulations employing harmonic con- 
straint potentials. The corresponding error is analogous 
to the convolution errors discussed in this paper. Note 
that the parameter a can be converted to a force con- 
stant for a harmonic potential via k = 2fc^T / 'a 2 , where 
k is the force constant. Errors resulting from finite k 
in TI and US computations can be identified as result- 
ing from a convolution between the true density of states 
and a known Gaussian function. Typically, the harmonic 
constraints are tight enough for A a to be a good ap- 
proximation of A but any persisting bias can, at least in 
principle, be removed by deconvolution as shown below. 



B. Interest of the mollified free energy 

In this work, we use A a to define an adaptive bias. 
The first interest of this approach is that the gradient of 
A a is much easier to compute than the gradient of A. 
Indeed, the laster reads (see References [10l - [l2l ]) 



( 4 ) 

where (•}£* denotes a canonical average for a fixed value 
of the reaction coordinate, and G is the Gram ma- 
trix. The latter matrix is defined as G — J J 1 with 
Jij = d£i/ 'dxj {xi are the Cartesian coordinates on which 
the reaction coordinates are defined). The computation 
of the free energy gradient therefore requires the compu- 
tation of second derivatives of the reaction coordinate, 
which is cumbersome in many cases. The gradient of the 
mollified free energy has a much simpler expression: 



dA a (C) 

da 



= -k B T JX 



dq6 a ^{x)-C)e~ 0V{x) dx 



6 a (^)-e)e- mx) dx 



-, (5) 
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where j is a reaction coordinate index and 

J QL 

In particular, no derivative of the reaction coordinates 
are required. 

Another interest of the mollified free energy lies in the 
nonlocality of 5 a , which allows a single observation of £ 
to contribute to A a for a range of values £* , leading to a 
faster convergence. The question is then whether there is 
a range of a for which: (i) a is sufficiently large so that 
A a could be estimated with fewer samples than what 
would be required to compute A and (ii) a is sufficiently 
small, so that A a is close enough to A to efficiently bias 
the dynamics. We show in Section lTV Cl that a large range 
of a satisfies these two conditions on a paradigmatic test 
case. 



A new ABP method 



1. Construction of the method 

To compute approximations of (j3|) and ([5]) as time av- 
erages along a trajectory xt driven by the potential func- 
tion V(x), we first assume that x t is ergodic with respect 
to the canonical ensemble. We may take Xt as a solution 
to the Langevin equation driven by the potential V, for 
example. Using trajectory averages, ([3]) can be obtained 
as the following longtime limit: 

e -PA a (e,t) = z -i L + J* 6a{ ^ Xs) „ r) ds ^ > (6) 

where the normalization constant Z t is 

Zf= J n ( 1 + l ~ D ds^j d£ = |fi| + t. 



The normalization constant ensures that 



d£ = I 



at all times t > 0. Notice that we implicitely assumed 
that the reaction coordinate has values in a finite space 
il. This is indeed the case when angles are considered. 
For unbounded reaction coordinates, it is always possible 
to restrict the sampling to important values of £(x). In 
practice, the range of the reaction coordinate needs to be 
truncated anyway. 
From ([S]), we obtain 



(8) 



J 1+ / 6 a (Z(x s )-Ods 

Jo 

which, in the longtime limit converges to ([5]) 



Now, a simple ergodic average such as or ([5J can 
of course not be used in practice since the dynamics at 
hand are usually metastable for complex systems, and 
the convergence of the time averages ^ and (J5J is very 
slow. We therefore need to bias the dynamics in order to 
remove the metastability. 

In what follows, we will consider a trajectory xt ob- 
tained from the equations of motion with the biased po- 
tential V+Vb- The idea behind adaptive method is to use 
the opposite of some current approximation of the free 
energy as a biasing potential, and to update the estimate 
as time goes on, in a way such that the bias eventually 
converges to the correct free energy. Here, we consider 
an adaptive biasing potential method, defined through 
the following update of the biasing potential V&: 



(9) 



where the renormalized current approximation of the 
mollified free energy e - ^ Aj4a( ^'*) is 





-PA a (£,t) 




max 









The parameter c in © is an important quantity in our 
method, which allows to tune the convergence rate of the 
method. We discuss its choice in Section III C 31 With 
these definitions, V& = — A a up to an additive constant 
which is chosen such that max[Vb] = c. Similarly, AA a = 
A a , again, up to an additive constant which is such that 
min[AA Q ] = 

Departing from standard ABP/ABF frameworks we 
use ideas from importance sampling to write (JB]) and ((5J) 
as time averages over biased trajectories 

-PA a {C,t) _ y-l 



V (l + J 5 a (^x s )-ne^ x ^ds\ 



( 10 ) 

where Z t is still a normalization constant ensuring ([7)1. 
and 



(7) dAa(C,t) 



da 



= -k R T- 



1+ / 5 a (C(x s )-C)e^ {Xs) ' s) ds 
Jo 

(11) 

The ABP method we discuss here is based on the bias- 
ing potential ([9]), updated with the current estimate of 
the free energy (flU)) . New configurations are obtained 
by integrating in time the biased equations of motion us- 
ing the simple estimate (jTTJ) for the biasing force. The 
convergence of this method is discussed in Section III C 31 
In fact, Eq. ([TO]) is a way to evaluate the convolution in 
Eq. © at each point £* using a biased trajectory. This 
gives us a precise understanding of how using finite a 
introduces error in the estimate A a and how to remove 
that error. This is a strength of our method which makes 
it unique. If we try to draw analogy with metadynam- 
ics, the framework of (|I0[) would imply the continuous 
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deposition of the Gaussians S a at each point £(xt) along 
the trajectory. Notice that in this analogy the Gaussians 
would be added to the density of states rather than to 
the bias potential, precluding us from going any further 
with the analogy. 



2. Time-discretization 

Let us briefly discuss the time-discretization of the 
method based on (j!?l- (fTI^ - ([TT]) . Assume that we have 
a suitable discretization where time is broken into parts 
of duration At so t — nAt and XiAt is written Xi. The 
biasing potential is now updated as 



e -P&A a {£,n) e f3c 



(12) 



where e ^ Ayl ^«<™) = e -M»K,f.)/ maXf ,[e-^«(5*,u)] 
and (fTUj) and ([Til are respectively replaced by 



(13) 



= z- i [i+Y d Sa(^i)-e)< 



i=0 



and 



dA a (C,n + l) 



da 



(14) 



£%*a«(iTi)-$ 



-kpiT- 



At t = we have exp[— j3A a (£, 0)] = 1/Zq. Let us empha- 
size again that the trajectory Xi is generated from biased 
equation of motion associated with the biased potential 
V + V b . 

The implementation only requires storing the current 
value of the numerator and denominator of Eq. (TT4")) at 
the points £*. In particular, Z n is never needed in prac- 
tice, (see Appendix [AJ The biasing force — Wb, needed 
for instance to integrate the Langevin dynamics, is ob- 
tained through Eq. (|T4"|). 



tant constant), so that (| 10[) leads to 

i+ [ ^({(i.l-D^Wi 

= lim 



fi+ / ^(^)-e')e m( * ,),8) J^ 

o V Jo / 

^(e(x)-r)e^ (yw+c) ^ 



,v 



« a K(a;)-Oe-^ (a)+0) d!Bde' 

The fact that, if a limit exists, then it is the correct one, 
is an important consistency check of the method. How- 
ever, we were not able to prove that the biasing potential 
indeed converges (this issue arises in all ABP methods 
while such an alalysis can rigorously be done for some 
ABF methods [HI]). 

Let us now look more carefully at the first iterations 
of the algorithm, in order to understand the role of the 
constant c in © or (fl"2")) . We base our considerations 
on the numerical discretization (|13[) to simplify the argu- 
ment. First, recall that the constant c does not change 
the longtime limit of the algorithm. However, it helps 
accelerating the convergence during the initial transient 
regime. The first iteration of (Q~3|) indeed shows that 



l + 6 a (0)ef>* 



When c is such that e^ c is small, 1) is raised by a 

small amount and the gradient of Vb encourages trajec- 
tories to move away from £* to some small extent. By 
increasing the value of c, we obtain a bias potential that 
pushes trajectories away from £* more strongly, hence in- 
creasing the efficiency of the bias potential, in particular 
at the early stages of the process. We therefore conclude 
that the value of c should be as large as possible while 
maintaining numerical stability. Not all ABP methods 
update their biases according to this rule, see the compar- 
ison between our approach and the standard Self-healing 
umbrella sampling algorithm in Section [III Al 



III. COMPARISON WITH OTHER METHODS 



3. Convergence and consistency 

It can be checked that, if the biasing potential 14 
converges in the long-time limit, then it converges to 
—A a up to an additive constant. Indeed, denoting by 
A a (0 = lh ri n _ >+00 A a (t;, t), the trajectory Xj is sampled 
according to the limiting canonical measure associated 
with the potential V — A a + C (where C is an unimpor- 



A. Self-healing Umbrella sampling 

Self-healing umbrella sampling [5] (SHUS) can be seen 
as a special case of the method presented here. SHUS 
can be written in terms of Eq. (1121) using the following 
time-dependent constant: 



e /3c(n) = max 



-/3A„(r,n) 
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With this choice we have, Vb = —A a and J n e^^'^ d£ = 
1. This choice for c(n) was suboptimal since the analy- 
sis of section III C 31 shows that the value of c should be 
as large as possible. Notice also that when the reac- 
tion coordinate space is discretized into a finite number 
of bins, the normalization condition ([7]) should be re- 
stated as a sum over bin indexes and the maximal value 
of exp[— f3A a (£, n)] is therefore less than one. This cor- 
responds to a negative value of c(n). We checked for 
the testcase considered in Section IIVI that our method 
outperforms SHUS for precisely this reason. 

B. Adaptive biasing force 

We compare numerically our approach to two ABF for- 
mulations in Section TlV Al ABF is a good reference for 
comparison because there are no model parameters to 
choose. Errors arise only through time and reaction co- 
ordinate discretization. Two exact formulations of the 
free energy gradient are (j4j above, and 

**•>--<=(*§)>«.• (15) 

where MjT 1 = JMJ 1 with M the mass matrix and J de- 
fined in Q (see reference [1] for this second expression). 
We point out that in practice F(£) is approximated by 
a trajectory average t) which is then used to bias 
the dynamics. For further details on the expressions Q 
and (|15[) or their numerical implementation, we refer the 
reader to the cited works. 

With ABF, one must address constructing the free 
energy from an estimation of its gradient, the calcu- 
lated field F. While there are specific solutions to this 
problem[H [3, HH we employ a standard variational for- 
mulation. We recast this question as an optimization 
problem where the objective function 

I(u)= [ ||F(0-V € u|| a d£ (16) 

is to be minimized. The corresponding Euler-Lagrange 
equation is 

A c u(0 = V £ • F, (17) 

which is just Poisson's equation, to be supplemented with 
appropriate boundary conditions (depending on the do- 
main f2). The solution u(£) is the best representation 
of the free energy given the vector field F(£). This is 
solved via finite difference in the present work, but finite 
elements (or any Galerkin method) could be used as well. 

C. Metadynamics 

Because we have developed a method within the adap- 
tive bias potential paradigm, we also make a comparison 



to well-tempered metadynamics [lfj. In this formulation 
of metadynamics the bias potential in one dimension is 
given by 

VT eta (£,T) = £ G((-&,h(£,t'),w), (18) 

t'<T 

where the functions G(X, H, W) are Gaussians of width 
W and height H, centered on X. We write y & meta 
to indicate that this is the bias potential generated 
by metadynamics. The Gaussian height in well- 
tempered metadynamics is both dependent on time 
and position along the reaction coordinate /i(£, t) — 
wexp[-V & meta (^,i)/fc B AT]r G . For details of this version 
of metadynamics we refer the reader to reference • We 
compare to this particular formulation because it requires 
less interaction with the user and a choice of parameter 
values is given in the cited reference. 



IV. NUMERICAL EXAMPLES 
A. Simulation details and results 

Alanine dipeptide is a familiar system for benchmark- 
ing sampling methods]!, 0, H3jla |. Here, we employ 
AMBER with a half femtosecond timestep, no con- 
straints, solvent effects are modeled with the generalized 
Born model and we use the ff94 parameterization. The 
temperature was maintained at T = 300K with Langevin 
dynamics where the collision frequency is 1 ps _1 . We 
select the common backbone dihedral angles (£1,62) = 
($, \&) as reaction coordinates. 

When discretizing the reaction coordinate, it is com- 
mon to use a small bin size to be sure that the free energy 
is correctly captured. Here, we use 300 bins of width 1.2 
degrees. We will also consider a bin width of 3.6 degrees 
for Eq. (fT5l) to examine the influence of bin size on ABF. 
In practice, for Eqs. (|13|) and (Q~4|), the current configura- 
tion along a trajectory may contribute only to an m by 
m grid centered around (£i(xt), £2(2^))) which amounts 
to truncating the range of the Gaussian function S a . The 
number of bins m were chosen so that 6' a (£ — £*) is negli- 
gible for £* outside this box. For example, when a = 5° 
we use m = 20. In practice we neglect the normaliza- 
tion Z n as well as the normalization of 5 a . We give a 
schematic algorithm in appendix [A} 

In simulations with equations (flU)) . (g| and (fT5|) . we use 
a "ramp function" R(Ni t k) — min[l, Ni k/ No] to scale the 
biasing force (see for instance reference [8j), where Ni t k is 
the population in bin (/, k) and the parameter A*^ = 10 
was optimized for equations (j4|) and (|15jl . The ramp 
function scales the biasing force so that the initially noisy 
observations of the force do not induce non equilibrium 
effects. The biasing force for the method presented here is 
given by Eq. (fTi)) . The biasing force for the ABF methods 
are given in equations ^ and (|15p . The biasing forces 
(and biasing potential) are updated at each timestep. 
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To study sampling efficiency we use the average differ- 
ence 



n n 

d w = ^EEK*( fc >o-^(M,*)| ( 19 ) 
fc=i 1=1 

between the estimated free energy and a reference to be 
defined below, n is the number of bins in each coordinate, 
k and I are bin indices. For Eq. (4} and (T"5} A is the 
solution of Eq. (TT} . In out method, A is either the left- 
hand side of Eq. (TB"} , A a , or its deconvoluted version 
A^ cnvl . Finally for Eq. (JTHJ), A = — (T + Ar)V h mcta /Ar. 
The reported results for d(i) are found by using only a 
single trajectory with each method. We do not report 
the results obtained with SHUS since the convergence 
was found to happen much slower than for cases where 
c> 0. 

We use the Richardson-Lucy algorithm [2(| [5l[ to 
deconvolute A a because of its simplicity but another 
method of deconvolution could be used, in particular if 
S a is not defined as a Gaussian. This algorithm is de- 
noted by "RL" throughout. The RL algorithm uses the 
following iterative procedure: 



(20) 

where /o(£) = exp[— /3A Q (£)], which is given by equa- 
tion (p~3|) . To begin the algorithm, S a and /$ must be 
normalized. The fixed-point iteration (|2U|) suggests that 
/«(0 ~~ * ex P[ — as n — > +oo. We use 10 iterations 

in the reported results. 

The reference free energy was computed by reweighting 
a long biased trajectory (120 ns) as 

A Ici (k,l) = ~k B T\n \J2 S [Zk -ti(xi)]x 
5^-6(^)]exp[/3H(C(^),r)] 



where Vh(£, r) was constructed from 4 ns of simulation 
with the mollified DOS method. The bias was not up- 
dated during construction of the reference free energy. 
This produces a result free from errors associated with 
the choice of a. The reference profile A re { is shown in 
figure QJA) and in figure HJB) we show A a (£,t) at 1 ns 
of sampling with a — 5°. The average difference d{t) 
is shown in figure [2] for A a with different values of a. 
To show how the zero of energy of the bias potential 
controls the speed of convergence, in figure [H we plot 
Eq. jT3} with c = in (12} and we set c = 15k B T for 
the remaining simulations. In figure [3] we show d(t) for 
Eqs. (gj) and (15} (ABF methods). Results for Eq. (15} 
(well-tempered metadynamics) are shown in figure |U[D). 



B. Efficiency of the results as a function of a 

For small a the nonlocality of the formulation dis- 
appears and in figure |5] we see slow convergence for 
a = 0.8°. For intermediate values of a, nonlocality al- 
lows the bias potential to equilibrate much faster. For 
a = 2° and a — 5°, A a is a good approximation of A, 
d(t) falls well under 1 kcal/mol and we observe high effi- 
ciency. With the value a = 10°, d(t) plateaus at roughly 
1 kcal/mol; a is now too large for A a to be a good ap- 
proximation of A. After applying the RL deconvolution 
to A a= iQo , d(t) drops to match the accuracy obtained 
with a — 2° or a = 5°. The correspondence between A a 
and A has deteriorated but not enough to decelerate the 
sampling: A a= io° is still a good biasing potential and A 
can be recovered with deconvolution even at very short 
times. 

For large a Eq. (IT41 approaches zero, leaving only a 
small biasing force to accelerate the dynamics. To assess 
whether a = 20° is so large as to slow down the sampling, 
we apply the RL deconvolution. The results in figure [2] 
demonstrate that A can be recovered to high accuracy 
for a — 20° at long times but that sampling efficiency is 
affected. 

In figure[3]we show d(t) for Eqs. (4j and (IT51) with a bin 
size of 1.2° and also for Eq. (15} with a bin width of 3.6°. 
If we compare the time to reach d(t) = 1 kcal/mol, sim- 
ulation with Eq. (13} is roughly three to ten times faster 
than Eqs. g} and (15} for 2 < a < 20 at the bins size of 
1.2°. For the larger bin size 3.6°, ABF sampling speed 
becomes competitive with the mollified density of states 
approach but it is impossible to remove the error. The 
3.6° bin width coincides with the Gaussian half-width of 
5 a when a = 2°. A larger bin size can enhance sam- 
pling speed for ABF but at a cost in accuracy. Note that 
a = 20° corresponds to a S a with a half width that spans 
33.3° in one dimension. This is a very large effective bin 
width for the accuracy of the results; A similar bin size 
with Eqs. d} or (T5} would produce large, irreparable 
errors. 

In figure U we show results for the metadynamics sim- 
ulations. We use the values AT = 1800 K, u = 0.24 
cal mol -1 fs -1 and tq = 120 fs, as suggested in refer- 
ence |16| . We could not improve the results by choosing 
different parameters. In panels (A) and (B) of figure @] 
we show the absolute difference between the computed 
A a= io and the reference A with and without deconvo- 
lution, respectively. Clearly, the bulk of error is due to 
the missrepresentation of the very negatively curved re- 
gions of the free energy and the ability to deconvolute 
A a drastically reduces this error. In panel (C) we show 
the absolute difference between the free energy computed 
via equation (15} and the reference. We see again that 
the error is concentrated in the regions of large negative 
curvature but there is not simple and obvious way to re- 
duce these errors with some post-process. Panel (D) con- 
firms that the metadynamics promotes extremely rapid 
sampling but that the long time accuracy, especially in 
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d(t) (kcal/mol) 



d(t) (kcal/mol) 




FIG. 2: Error (|19p as a function of time for the method pre- 
sented in this paper, for various a with and without decon- 
volution. Unless otherwise stated, c = 15/cbT. The a — 5, 
c = simulation demonstrates slow convergence due to a sub- 
optimal choice of c, as described in the text. In the inset we 
show the last 4 ns of the a = 5° results. 



strongly curved regions, is limited. 

The results summarized in figures[2l[3]and[4]imply that 
a wide range of values 2 < a < 20 lead to good efficiency. 
The ability to use the simple deconvolution algorithm is 
a clear advantage of the method. 



C. Choosing a a priori 

We now discuss how to a priori choose a based on 
some rough error estimates. Taking £ as a scalar, we 




a-5 
RL a=5 
Eq4 
Eq 15 
Eq15(100 bins) 



1.5 



2 2.5 
t(ns) 



3.5 



FIG. 3: Error (|19|) as a function of time for a = 5° and a = 
10° with the method proposed in this paper, and comparison 
with ABF results obtained from Eqs. © and (fT5)l . 



may expand e ) as a Taylor series. Eq. ^ yields 



A'(£*)y A"(C) 



(21) 



where we keep terms up to the second moment of S a . 
Assuming that A(£) is harmonic near the minimum £ = 
q, the curvature can be estimated as A"(q) — ksT/a 2 
where a 2 is the variance of the reaction coordinate at 
temperature T. From Eq. (|2Tj) . 

exp[-P(A a (q)-A( q ))]~l-^. 

While the higher order terms and the regions where 
A' are certainly important to the total error, this 
motivates defining a as a function of a if little is known 



8 



|A a=10 -A| (kcal/mol) |A a=1 dcnvl -A| 




-150 -100 -50 50 100 150 -150 -100 -50 50 100 150 




FIG. 4: In panel (A) we show the absolute difference between the computed A a =io and the reference A. Most of the error is, 
as expected, due to the regions of large curvature. In panel (B) the absolute difference between the deconvoluted free energy 
AtTio and the reference A is shown. In panel (C) we show the absolute difference between A meta = -(T + AT)V b meta / AT 
and the reference A. The free energy estimates in panels A-C were taken at the end of a 4 ns trajectory. In panel (D) (|19[) is 
shown for the well-tempered metadynamics method (|18[) for comparison with results from the mollified DOS method. 



about the free energy — we can always calculate a in the 
initial state. 

We calculate the variance of the reaction coordinates 
to be about a 2 = 340 degrees squared for both $ and 
vF In terms of the values of a discussed above, this 
implies a/9 < a < a/2 as a good range for fixing a from 
calculation of a. Of course, different a's may also be used 
for different coordinates. 



deconvolution can be applied at the end of a simulation to 
remove all of the errors associated with the choice of the 
model parameter a — a unique feature and strength of 
this approach. This is limited only by the extent of sam- 
pling and the spacial discretization. This scheme easily 
accommodates the computation of the free energy sur- 
face and free energy gradient in several dimensions. We 
also suggest a simple means of a priori specifying a and 
c that should be quite general in applicability. 



V. CONCLUSION 



In conclusion, we have developed and tested an effi- 
cient ABP scheme. The nonlocality of 6 a leads to a bias 
potential and a bias force that equilibrate rapidly. Shift- 
ing the zero of energy on the bias potential was shown to 
result in efficient importance sampling. The parameter 
c has influence on only the efficiency of the importance 
sampling but not on the limiting error of A a . Because 
the bias potential is related to a convoluted free energy, 
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Appendix A: A Schematic Algorithm 

To help illustrate the simplicity of implementing equa- 
tion (7) from the text 



dA a (£*,n + l) = i=g 



X;%*«(e(»i)-r)e^ (€(B<)ii) 



'l.TO!£(K),i) 



i=0 



( A1 ) 

for a 2 dimensional computation, we give a schematic 
algorithm here. We first define some array names. Let 
the array named "pop(fc,l)" store the population at the 
(k, I) grid point (this is just the denominator of Eq. (|A1[) 
above), where k corresponds to the bin index of 
and I corresponds to the bin index of £2(2^)- Let the 
array named "dpop(j, k, /)" hold the derivative of the 
population along the £j=i,2 direction at the point (fc, /). 
The array "dpop" is simply the numerator of Eq. (|Alj) 
above. We use "dA(k)" to store the gradient of the free 



energy at the present point (k,l). We assume that a 
has been calculated and c has been specified. We let 
5 a (t;) = e~'^ I a , which amounts to ignoring the nor- 
malization of the Gaussian functions. Lastly, we denote 
the trajectory in phase space as Xi, F{n') is the force 
along the n' th degree of freedom, d/dn! is the derivative 
with respect to the n' degree of freedom and we use V (x) 
for the potential energy. 

First we initialize the arrays. 

t = 0, pop(fc,Z) = l V k, I and 
dpop(j, k, I) = V k, l,j and M = 1, 

where M = max^; [pop(fe, I)]. Each time the molecu- 
lar dynamics forces are computed we must also compute 
the current biasing information. Notice that we define 
everything in terms of the "pop" and "dpop" arrays so 
that no array is needed for the bias potential and that 
M = max / t i ;[pop(fc, /)] can be updated without looping 
over the full reaction coordinate domain. 



! evaluate free energy gradient at (k,l) for j = 1,2 
dA(j) = dpop(j,fc, Z)/pop(fc, I) 

! add bias forces to the existing forces and use a 

! ''Ramp function'' R as described in the text 
R = min(l,pop(yfc,0/10) 
F(n') = F(n') + R £* =1 dA(j) dfc/dn' 

! evaluate the weighting factor W for updating ''pop'' and ''dpop'' 
W = exp[£ V b ] = exp[/3 c]pop(fc,0/M 

! update ''pop'' and ''dpop'' on an m by m grid 
loop k! — k — m/2, k + m/2 
loop I' = 1- m/2, 1 + m/2 

P op(fc',n= pop(*M0 + - £i* fe ')<U6M) - Q,i>)W 

if pop(k',l') > M then M = pop(fc', V) 
loop j=l,2 

dpop(i,fc',0= dpopO',fc',n+%[^(eiM) - Q,v)6 a fa(xJ) - Q tl ,)]W 



We have defined k' and V so that "pop" and "dpop" are 
updated on an m by m grid as discussed in the text. The 
treatment of (k', I') should reflect whether the domain 
is assumed to be periodic or not. The approximate free 
energy A a is recovered (up to an additive constant) with 
A a = k B Tln[pop{k,l)/M]. 

The dynamics will now evolve in the presence of the 
biasing force dA(j), while the arrays "pop" and "dpop" 



hold unbiased estimates of the population and the deriva- 
tives of the population. Notice that the free energy gradi- 
ent is reduced to a simple ratio and the only difficulty lies 
in the careful treatment of the loops over the grid points 
k' and V . The often mathematically complex computa- 
tion of the free energy and free energy gradient is reduced 
to simple bookkeeping. 
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